Hofstadter butterflies of bilayer graphene 
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We calculate the electronic spectrum of bilayer graphene in perpendicular magnetic fields non- 
pert urbatively. To accomodate arbitrary displacements between the two layers, we apply a periodic 
gauge based on singular flux vortices of phase 2ty. The resulting Hofstadter- like butterfly plots show 
a reduced symmetry, depending on the relative position of the two layers against each other. The 
split of the zero-energy relativistic Landau level differs by one order of magnitude between Bernal 
and non-Bernal stacking. 
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After the theoretical prediction of the peculiar electronic 
properties of graphene in 1947 by Wallace 1 and the sub- 
sequent studies of its magnetic spectrum^ it took half 
a century until single layers of graphene could be iso- 
lated in experiment^ and the novel mesoscopic properties 
of these two-dimensional (2D) Dirac-like electronic sys- 
tems, e.g., their anomalous quantum Hall effect, could 
be measured.— Inspired by this experimental success, 
graphene has become the focus of numerous theoreti- 
cal works i 8 ! 9 ! 10 ! 11 ! 12 For bilayer s of graphene, an addi- 
tional degeneracy of the Landau levels and a Berry phase 
of 2tt were predicted to lead to an anomalous quantum 
Hall effect, different from either the regular massive elec- 
trons or the special Dirac-type electrons of single-layer 
graphene^ which was confirmed in experiment shortly 
afterwards^ and used for the characterization of bilayer 
samples J£ 

The low-energy electronic structure of a single layer of 
graphene is well described by a linearization near the cor- 
ner points of the hexagonal Brillouin zone (K points), 
resulting in an effective Hamiltonian formally equivalent 
to that of massless Dirac particles in two dimensions^ 
A related Hamiltonian can be constructed featuring a su- 
persymmetric structure which can be exploited to derive 
the electronic spectrum in the presence of an external 
magnetic field J£ The level at zero energy, characteristic 
for any supersymmetric system, maps directly to a spe- 
cial half- filled Landau level fixed at the Fermi energy £p, 
henceforth called the supersymmetric Landau level (SU- 
SYLL). 

In this Rapid Communication, we use the nonperturba- 
tive method pioneered in 1933 by Peierlsi^ for the im- 
plementation of a magnetic field in a model, which led 
Hofstadter, in 1976, to the discovery of the fractal spec- 
trum of 2D lattice electrons in a magnetic field. 19 Since 
its discovery, various aspects of the so-called Hofstadter 
butterfly have been studied , 2Q i 21 particularly in relation 
to graphenelike honeycomb structures i 12 i 22 i 23 Featuring 
a large variety of topologies, all these systems have in 
common that the atoms inside the unit cell are located 
on discrete coordinates. All closed loops have commen- 
surate areas, and the atomic network is regular enough 
that the magnetic phases of all links can be determined 
individually without the need of a continuously defined 
gauge field. For bilayer graphene, such a direct scheme 
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FIG. 1: (Color online) Hofstadter butterfly of a bilayer 
graphene in the Bernal stacking configuration. The band 
structure at zero magnetic field is rotationally symmetric in 
good approximation for an area around the K point and shows 
a split into four massive bands, with the two middle ones 
touching at Ef. The density of states (DOS) of a finite- width 
ribbon [a pair of (200,0) zigzag ribbons] in the same configu- 
ration shows the SUSYLL emerging at finite magnetic field. 
The split of the SUSYLL (discussed below) is not visible due 
to the limited resolution of the plot. 



for implementing a magnetic field is possible only for 
highly symmetric configurations like Bernal stacking ! 13 i 24 
To handle more general configurations, such as contin- 
uous displacements between the layers, it is in general 
unavoidable to choose a continuously defined gauge that 
fixes the phase for arbitrarily placed atoms. The difficulty 
that arises can be seen immediately: For any gauge field 
that is periodic in two dimensions, the magnetic phase 
of a closed loop around a single unit cell must cancel out 
exactly, corresponding to a vanishing total magnetic flux. 
Conversely, this means that any gauge field that results 
in a nonzero homogeneous magnetic field will invariably 
break the periodicity of the underlying system. 
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A possible way to bypass this problem is based on defining 
a magnetic flux vortex, here oriented in the z direction 
and located in (xo?2/o)> a ^ 2 ^ 26 

B(x,y,z) = $ 5(x-x )5(y-y )e z: 

where <£>o = h/e is the flux quantum. Physically, such a 
vortex is equivalent to a vanishing magnetic field, since 
it leaves the phase of any possible closed path unchanged 
modulo 2tt. One possible gauge field resulting in such a 
single flux vortex can be written as 
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Finding a periodic gauge follows straightforwardly. To 
the homogeneous magnetic field, we add a periodic ar- 
ray of flux vortices with a density such that the average 
magnetic field is exactly zero. For the resulting field, 
which is physically equivalent to the original, it is now 
possible to find a gauge field with the same periodicity 
as the array of vortices. If the underlying system is pe- 
riodic and the array of flux vortices has commensurate 
periodicity, there exists a supercell where the magnetic 
Hamiltonian is periodic. One possible periodic gauge that 
is especially advantageous for numerical implementation 
consists in a two-dimensional periodic system with lat- 
tice vectors a x and a y . The reciprocal lattice vectors 
(scaled by 2tt) are di such that ai • dj = 5ij . The mag- 
netic field is B = £$>o (a x x d y ) with I integer. The 
usual linear — but aperiodic — gauge for this field would 
be Aii n (r) = (r • d x )d y . A periodic gauge can now 
be defined as: 

A(r) = e$olr-a x }(a y -6(lr-a y })a x ) 

where [•] denotes the fractional part of a real number. 
To make sure that the phase of every link between two 
atoms is well defined, the gauge field is displaced by an 
infinitesimal amount such that every atom is either left 
or right of the divergent line. 

The Hamiltonian without magnetic field — based on a 
tight-binding parametrization originally used for multi- 
walled carbon nanotube o 23 i 27 — consists of a contribution 
for nearest neighbors within a layer and one for pairs 
of atoms located on different sheets 
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In absence of a magnetic field, the intralayer hopping 
is fixed to 7^ tra = 70 = 2.66 eV, while the interlay er 
hopping depends on the distance only, 



(3 exp 



with f3 = 70/8, a = 3.34 A, and S = 0.45 A. A cut- 
off is chosen as r cuto ff = a + 56. Following the Peierls 
substitution^ the magnetic field B is now implemented 
by multiplying a magnetic phase factor to each link be- 
tween two atoms i and j : 

litj (B) = 7iJ (B = 0) exp / 3 A B (r) • dr 
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FIG. 2: (Color online) Hofstadter butterfly of a bilayer 
graphene in two differently shifted configurations. Top panel: 
A A stacking (two layers exactly aligned). The band structure 
for this highly symmetric stacking (same rotational symmetry 
as for Bernal stacking in Fig. [T]) shows the single-layer cone 
simply split up in energy. Bottom panel: Intermediate posi- 
tion between Bernal and AA stacking. The rotational sym- 
metry is broken and the bands split into two cones at different 
offsets from the K point and different energies. The straight 
lines overlaid at the energy minimum and maximum are the 
regular Landau levels of the massive bands. Near Ef, one can 
make out the parabolic traces of the relativistic Landau levels 
and the horizontal lines of the SUSYLLs (see text). Insets at 
the lower right of each panel: DOS of a finite-width ribbon 
shows the corresponding behavior in each case. 
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FIG. 3: (Color online) Evolution of the split of the supersym- 
metric Landau level as a function of the displacement between 
the two graphene layers. Top panel: Magnitude of the split for 
displacements in two directions. The light spots correspond to 
Bernal stacking where the level is nearly degenerate. Bottom 
panel: Same data along a cut at Sy — 0. The small remain- 
ing split at the Bernal stacking configuration originates in the 
long-range interlayer hoppings contained in the parametriza- 
tion. The small discontinuities are caused by the cutoff r cu t- 
The calculation here was done at <I> = <J>o/256, but proved to 
be independent of the magnetic field for values up to ~ 0.05$o • 

where the integral is computed on a straight line between 
the atomic positions and rj. 

For the bilayer graphene, we arrive thus at a periodic 
Hamiltonian with a two-dimensional unit cell contain- 
ing four atoms and spanning the area of one hexago- 
nal graphene plaquette: A p i aque tte = (3>/3/2)d£ c , where 
dec = 1.42 A is the intralayer distance between neighbor- 
ing carbon atoms. The effect of a perpendicular magnetic 
field, measured in flux per plaquette <1> = A p i a q Ue tte^> 
can be calculated for commensurate values $ = (p/q) $0 
(p, q integers) by constructing a supercell of q unit 
cells. The corresponding Bloch Hamiltonian H (k) is a 
4g x 4g matrix that can be diagonalized for arbitrary val- 
ues of k in the two-dimensional Brillouin zone of area 
47r 2 /gApi aquet te- 

To obtain the butterfly plots as displayed in Figs. [1] and 
[2j we chose 0^p^g = 512, reducing the fraction p/q 
for efficiency. For each value of $ the density of states 
was calculated from a histogram over the spectral values 
for a random sampling of k over the Brillouin zone. The 
number of sampling points was chosen individually for 
different values of p to achieve convergence. In Figs. Q] 
and [2j the Hofstadter spectra of three differently aligned 



graphene bilayers are presented. The Bernal stacking 
(Fig. [1]) stands out, as it is the configuration of layers 
in natural graphite j 24 i 28 Alternative configurations like 
AA stacking were found in ah initio calculations to be en- 
ergetically unfavorable; 29 they can, however, be thought 
of as either mechanically shifted samples or sections of 
curved bilayers (e.g., sections of two shells in a large mul- 
tiwall carbon nanotube) where the alignment unavoid- 
ably varies over distance. Compared to the Hofstadter 
butterfly of a single sheet of graphene^ two asymmetries 
are visible in all three plots: The electron-hole symme- 
try (E <-> —E) is broken down by the interlayer coupling 
already at zero magnetic field: while the lowest-energy 
states of a single graphene layer have constant phase over 
all atoms and can couple efficiently into symmetric and 
antisymmetric hybrid states of the bilayer system, the 
states at high energies have alternating phases for neigh- 
boring atoms, so interlayer hybridization is prohibited 
by the second-nearest-neighbor interlayer coupling. For 
low magnetic fields, two sets of Landau levels can there- 
fore be observed at the bottom of the spectrum, indi- 
cating a split of the massive band of graphene at the Y 
point (-E^in = — ^7o, ttIq = 2h 2 /^qo1^ c ) into two bands 
at different energies and with different effective masses 
[^min » £min±l-l eV, m* ± « m* /(l T 2. 1/3/Tb) , indepen- 
dent of the relative shift of the two layers; see the straight 
lines overlaid in the bottom panel of Fig. [2]. At the top 
of the spectrum, where the split is prohibited, only one 
degenerate set of Landau levels appears, as in single-layer 
graphene. The original periodic symmetry along the B- 
field axis at one flux quantum per graphene plaquette is 
broken down due to the smaller areas formed by inter- 
layer loops. The breaking of this symmetry is compara- 
bly small in the ^4^4 stacking configuration (Fig. [21 top) 
where loops of the full plaquette area are dominant. In 
the two other configurations smaller loops are more domi- 
nant, so the periodicity is perturbed more severely. In the 
intermediate configuration (Fig. [21 bottom), the fractal 
patterns appear slightly smeared out for high magnetic 
fields, due to the reduced symmetry of the system. 

The right insets of Figs. [1] and [2] display the spectra 
of (200,0) bilayer graphene nanoribbons^U each in a 
corresponding configuration, obtained by a method de- 
scribed before 23 that allows handling of continuous mag- 
netic fields For low magnetic fields, these spectra are 
strongly influenced by finite-size effects*^ Only for mag- 
netic fields larger than B* w 4<£>o/<i 2 , which for a rib- 
bon of width d = 50 nm relates to ~ 7 T, do the spec- 
tra of two-dimensional bilayer graphene begin to emerge. 
Prominent in all three insets are the dark, horizontal pairs 
of lines at the center, the super symmetric Landau levels. 
While these represent discrete levels in two-dimensional 
graphene sheets, they are broadened by the finite width 
of the ribbon to a peak of the same shape as in carbon 
nanotubes j 23 i 32 The mesoscopic character of these split 
SUSYLLs in dependence on the width W of the ribbon is 
captured by the functional form of the density of states: 

p(E,B,W) = f((E-E )W,BW 2 ) 

where Eq is the position of the maximum. 



Single-layer graphene is known to feature an anomalous 
supersymmetric Landau level at the Fermi energy^&ii 
Neglecting Zeeman splitting, this level is fourfold degen- 
erate (twice spin, twice valley) and half filled. For bilayer 
graphene in Bernal stacking (Fig. [T|) the SUSYLLs of the 
two layers have been shown to be protected by symme- 
try and to remain degenerate, giving in total an eightfold 
degeneracy^ In Fig. [21 this degeneracy can be observed 
to be lifted for displaced bilayer s, leading to a split of the 
SUSYLL into a bonding and an antibonding hybrid state 
in the two layers, each fourfold degenerate. The continu- 
ous evolution of the split for varying displacement of the 
two layers against each other is displayed in Fig. [3j The 
split reaches its maximum of AE ~ 0.3 eV for the ^4^4- 
stacking configuration and is minimal for Bernal stacking. 
For simpler tight-binding parametrizations that take into 
account only first- and second-nearest-neighbor interlayer 
hoppings, the degeneracy in the Bernal configuration is 
known to be exact. 13 Here, in contrast, this degeneracy 
is split by AE ~ 0.01 eV due to interlayer hoppings of 
a longer range, similar to the effect caused by second- 
nearest-neighbor interactions within one layer. 33 



In conclusion, we have developed a method that allows 
the nonperturbative implementation of a magnetic field 
in periodic systems with arbitrarily positioned atoms. A 
7r orbital parametrization for graphitic interlayer inter- 
actions with arbitrary displacements was then used to 
calculate the Hofstadter spectrum of bilayer graphene 
in various configurations, revealing common features like 
electron-hole symmetry breaking, and differences, espe- 
cially in the breaking of the magnetic-field periodicity. 
A close look at the supersymmetric Landau level at low 
fields near the Fermi energy revealed a breaking of the 
previously found symmetry, resulting in a split of the 
level, depending on the lateral displacement of the two 
graphene layers against each other. 
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